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We study the flow behind an array of equally spaced parallel cylinders. A system of Stuart- 
Landau equations with complex parameters is used to model the oscillating wakes. Our purpose 
is to identify the 6 scalar parameters which most accurately reproduce the experimental data of 
Chauve and Le Gal [Physica D 58, pp 407-413, (1992)]. To do so, we perform a computational 
search for the minimum of a distance J7\ We define {J as the sum-square difference of the data 
and amplitudes reconstructed using coupled equations. The search algorithm is made more 
efficient through the use of a partially analytical expression for the gradient Vj7\ Indeed VjT 
can be obtained by the integration of a dynamical system propagating backwards in time (a 
backpropagation equation for the Lagrange multipliers). Using the parameters computed via 
the backpropagation method, the coupled Stuart-Landau equations accurately predicted the 
experimental data from Chauve and Le Gal over a correlation time of the system. Our method 
turns out to be quite robust as evidenced by using noisy synthetic data obtained from integra- 
tions of the coupled Stuart-Landau equations. However, a difficulty remains with experimental 
data: in that case the several sets of identified parameters are shown to yield equivalent predic- 
tions. This is due to a strong discretization or "round-off" error arising from the digitalization 
of the video images in the experiment. This ambiguity in parameter identification has been 
reproduced with synthetic data subjected to the same kind of discretization. 
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1 Introduction 



When analyzing chaotic or turbulent dynamical systems one is often faced with 
the problem of choosing optimal models. For instance, there are several models 
of turbulence, such as the Reynolds stress, eddy viscosity or k-e models, that are 
frequently used by practitioners. One would obviously like to make an optimal 
choice among this wealth of models. There are two aspects to this: to find the 
general functional form of the model, and to find the parameters for a given 
form. In this paper we report our experience with such a parameter identification 
problem. 
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The experimental system consists in a grid of 16 evenly spaced cylinders 
placed transversally to a laminar incoming flow. For a single cylinder, when 
the Reynolds number Re = ^ — based on the upstream velocity U, cylinder 
diameter d and viscosity v — exceeds a critical value Re c ~ 46, the celebrated 
Benard-Von Karman vortex street sets in. In this regime the wake of the cylinder 
behaves like a nonlinear oscillator. When an array of cylinders is used as in 
ref. [1, 2, 3] one obtains instead a chain of oscillators. Such systems have been 
studied theoretically for some time now (for some recent references see [4]). There 
are many examples in which they display irregular oscillations just as seen in 
Figure la. The experiment of ref. [1, 2] thus provides a particularly simple 
hydrodynamical example of spatiotemporal chaos. 

The Hopf bifurcation of a single cylinder may be described [5] by a Stuart- 
Landau equation with complex parameters 

d t A = a A- l\A\ 2 A (1) 

where a = r + ku, / = l r + i/;, and A is a complex amplitude related to the 
observations of the cross-stream deviations of the wakes through its real part. 
This equation is a normal form for the Hopf bifurcation. Now in the experimental 
system, the individual wakes with amplitude Ai will undergo some degree of 
coupling with their neighbors. The simplest possible form for this coupling is a 
linear finite-difference operator [6] 

T d t Ai = a A, + g(A i+1 + A_! - 2Ai) - Z|A,fA t - (2) 

where g = g r + igi is a new complex parameter. Boundary and initial conditions 
will be discussed later. 

Although the three complex parameters <r, / and g may in principle be obtained 
from an analysis of the Navier-Stokes equations, in practice this is a formidable 
problem, requiring extensive numerical simulations even for a single cylinder. On 
the other hand, simply searching for the parameters that offer the best possible 
match to the data requires only small-scale computations with the methods we 
describe below. To give a more precise definition of our problem, we start with 
a set of 16 times series Vik = K'(^fc), where the index i indicates the cylinder and 
the tk are the measurement times. We connect the model solutions Ai(tk) and 
the real data Vik with the assumption 

v ik = MMh)) + W ik (3) 

where Wik is a decorrelated noise. Equation (3) expresses the ideal situation 
where the measurement process adds a Gaussian, independent error to each data 
value. In what follows however, we will discuss more complex measurement pro- 
cesses (see equation (29) below). Because of these doubts about the validity of 
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equality (3) it is useful to test our parameter identification procedure on a syn- 
thetic set of data generated directly from a numerical integration of equations 
(2) and (3). This use of synthetic data offers striking evidence of the role of ex- 
perimental artifacts such as large noise amplitudes, or other peculiarities of the 
measurement process. 

The selection of an appropriate parameter identification method is a still 
largely open question. While parameter identification, or more generally the tack- 
ling of inverse problems, is classical in various fields of science such as robotics or 
geophysical research, there is only one attempt we know of to deal with a spatially 
extended, chaotic dynamical system such as equation (2). Chauve and Le Gal 
[2] used the proper orthogonal decomposition (POD) to analyze data from the 
cylinder wake experiment. They were able to extract all the parameters except 
g r using the POD. With these parameter values and ad hoc values of g r they 
produced simulations that were roughly similar to the experimental data, but 
with incorrect phase relations. While the experimental data show that neigh- 
boring wakes are oscillating in opposition to each other, the POD yields wakes 
oscillating in phase with each other. 

In this paper we suggest another identification technique, which is based on 
an optimization procedure. The optimization consists of a search in the space of 
all solutions of equation (2) obtained when the parameters are varied. We define 
a distance J between the experimental data and a given solution. We then at- 
tempt to locate values of the parameters yielding the minimum J . This search is 
usually difficult, because the investigated space is large and local minima abound. 
Among the many methods that have been proposed, one of the simplest is to per- 
form a descent in the direction of the gradient which is estimated numerically. 
(Much more sophisticated methods may be found such as simulated annealing 
[7] or multigrid methods [8, 9]). We chose here to employ a more sophisticated 
descent method (the PLMA method described in the Appendix). The main dif- 
ference with the naive approach lies in the fact that the gradient is not estimated 
numerically. Instead it is obtained using a technique of backpropagation resulting 
in an adjoint dynamical system which propagates backwards in time. Although 
the the field of chaotic model reconstruction has been very active [10, 11, 12], the 
backpropagation technique is relatively new. We note that it has been used in 
ref. [13] for orbit reconstruction, but not for model reconstruction or, as we do 
here, for parameter estimation . 

In Section 2, we present the experimental setup and data. Section 3 introduces 
the theoretical model (basically equation (2)) that we take throughout this article 
as a working hypothesis. In Section 4 we describe the method used for parameter 
identification, including the definition of the cost function or distance J and 
we describe the backpropagation equations for the Lagrange multipliers. Finally 
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Section 5 contains our numerical results for real and synthetic data sets. 



2 Experimental Data 



Although the data set we analyze in this paper is different from the one used 
in [f, 2], there are no differences in experimental setup. The 16 cylinders are 
located in a hydrodynamic water channel. Each cylinder has a 2 mm diameter 
and is 200 mm long. The cylinder axes are 8 mm apart. Data were gathered 
above the onset of oscillations at Re = 80. Figure la gives a visualization of 
the 16 coupled wakes obtained by the emission of a dye streak at the back of 
each cylinder. Video images of the wakes are then digitized on 512x512 pixels. 
The line of pixels located 12 mm downstream of the cylinder's axes is recorded 
every T = 40 ms, yielding data similar to Fig. lb. These digitized lines are 
then processed to find the location of the dye streaks. In this first round of 
data processing, one wants the signal to pass near the maxima of intensity on 
the videoline. Pixels may have 256 intensity values. Thus there may be several 
relative maxima or adjacent pixels with the same maximal intensity. The value 
of the signal is then selected by a procedure that maximizes the smoothness of 
the signal. Details are given in [14]. The resulting deviation from the centerline 
position defines the data K'(^fc) plotted on figure 2a (with tk = kT). Some of the 
data used in this article will be available at the ftp site ftp.jussieu.fr in the 
directory jussieu/labos/lmm/Wakes. 



3 The Theoretical Model 

It is convenient to rewrite equation (2) in the form 

f = H(A;a) (4) 

where A(t) = (Ai(t))i = \ t N is a vector of N complex amplitudes, the parameter 
vector is noted a = (r, u, / r , g r , (/;), and H is the vector field defined by 

Hi(A; a) = (r + iu)Ai + g(A i+1 + A_ x - 2A 8 ) - Z|A,fA t - (5) 

Note that cylinders 1 and N pose a special problem since they lie at the ends 
of the chain. A simplifying hypothesis consists in defining 2 "virtual cylinders" 
outside the array of the N = 16 real cylinders for which a zero amplitude is 
imposed. Thus 

A = A N+1 = 0. (6) 
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The above equation plays the role of a boundary condition for our problem. 

To analyze the behavior of the system, let us however consider space-periodic 
solutions, neglecting for a while the effect of the boundary conditions. For r — 
4g r > and l r > 0, equation (4) has plane wave solutions of wave number q 

A t = Be l ^ t+q ^ (7) 

where 

fr + g r (2cosq-2) 
B= { 1 

and 

ft = u + g t (2 cos q- 2) -l t B 2 . (9) 

The wavenumber q should obey r + g r {2 cos q — 2) > . From these solutions it is 
easy to grasp the physical meaning of the various parameters. The distance to the 
instability threshold for in-phase (q = 0) oscillations is proportional to r, and the 
amplitude of the oscillations grows with it. The coupling between neighboring 
oscillators is a discretized version of the differential operator V 2 A that appears 
in the continuous-space Ginzburg-Landau equation. With positive g r this term 
tends to damp plane wave oscillations with a non uniform phase (q ^ mod27r). 
With g r < oscillations with a phase shift between cylinders are amplified, and 
the maximum amplification rate is achieved when q = ir. The parameters u, gi 
and U control the frequency of the oscillations. 

To perform computations, we need to discretize the continuous-time model 
(4). We hence divide the time window (t 0} t + jm&xT) into M — 1 time steps 
r. Experience shows that r must be much smaller than T, on the order of 
T/100, to reach accurate results. The discrete system is a simple forward Euler 
discretization 

A fc+1 = A fc + rH(A fc ;a) for k = 0, • • • , M - 1. (10) 
where = A(t + kr). 

Taking another point of view it is possible to chose a large time step r = T 
and to attempt to minimize the distance J with the resulting equation although 
it is an obviously poor approximation to equation (4). In this new point of view 
it is the discrete model rather than the continuous model which is postulated. 
This empirical approach is convenient because the discrete times correspond one 
to one to the data points, and that the computational work needed to integrate 
the model is scaled down by about 2 orders of magnitude 1 . 

1 It also allows for a simple type of least-squares fit of the parameters on which we will report 
elsewhere 
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4 Method 



To identify the parameters, we first need to define a distance J between exper- 
imental observations and model predictions. Let <fr(A ,t;a) be the flow inte- 
grated from equation (4) with parameters a, defined so that if A = A(t ), then 
A(t) = <fr(A ,t — t ;a) is the solution at time t. As in equation (3) the obser- 
vations are related to the real part of the amplitudes A(t). We let t°^ s = jT 
be the times at which observations are made, with T the sampling period of the 
experiment. We thus define the distance by 

imax , , 

J(A ,<i;t ,j max ) = J2 \\V(tj) ~ Re[<&(A ,tf s -t ; a )}\\ 2 (11) 

3 = 1 

where V(t) = (Vi(t))^z^ is the A-vector of time series of observations, jmaxT 
is the length of the series used for the definition of J7", and || . || 2 is the distance 
squared defined by ||X|| 2 = YliXf. In expression (II), the time t need not 
be the beginning of the full data set represented on Fig. 2a. Likewise, there is 
some leeway in the choice of the length jmaxT. Thus, we may define a distance 
for each time window (t 0} t + jm&xT). For each of these "training windows" the 
parameter identification is cast as an optimization problem, in which one searches 
for a global minimum of J with respect to a set of parameters. In the simplest 
version of this problem we express it as 

Find a op G R 6 such that J(A , a op ) < J(A , a) Va G R 6 

In more complex versions, we may add other parameters that we wish to vary, 
such as the initial conditions A . 

In order to solve this problem computationally, we use the PLMA algorithm 
proposed by Gill and Murray [15] and described in the Appendix. This method 
provides a minimum of J when its gradient V a j7" is available. While it is possible 
to find a gradient direction from several numerical estimates of J7", an efficient 
and accurate computation of the gradient is crucial since identification problems 
are generally ill-conditioned and lead to a large number of evaluations of the 
gradient. We thus construct an explicit analytical expression for the gradient 
V a j7\ This expression is derived from the problem defined in discretized time as 
in equation (10). This is a requirement which has been observed in practice to 
condition the success of the optimization algorithm. The distance thus reads 

M 

J^AcAv^AAf) = J2 e k\\V k -X fc || 2 (12) 

k=i 

where stands for the vector of experimental amplitudes at time kr, we de- 
composed the amplitudes into real and imaginary parts A = X + iY, and we 
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introduced an index e k which is equal to one when experimental data are avail- 
able — that is every T = 40ms — and zero in the opposite case. Finally we have 
used the notation J\ — instead of J — to express the fact that J\ is a function 
of the whole time series of amplitudes A , Ai, • • • , Am- When these amplitudes 
are related by the discrete model (10), we have J = J\. 

This completely defines the cost function in discrete time, but the dependence 
of Ji on a is implicit through its dependence on the variables A&. This difficulty 
can be handled using the Lagrange multipliers method. We hence recall a basic 
fact about constrained differentiation: 

Let Ji(a,x) be a function of R n X R m — > R. Consider the set of m constraints 
c(a } x) = where c is a function of R n X R m — > R m . Assume that the implicit 
function theorem applies so that these constraints define x as a function x(a) of a. 
Let J be the function of R n — > R defined by J(a) = J\(a, x(a)) . Call Lagrangian 
the function of R n X R m X R m — > R defined by 

C(a,x,p) = Ji(a^x) + c(a,x) ■ p, (13) 

where the components of p £ R m are called Lagrange multipliers. Under these 
conditions, if 

V P C = (14) 

and 

V X C = (15) 

then 

V a J = V a C. (16) 



We use this fact to compute explicitly the gradient of J . We introduce the 
multipliers = (Pik)i } Qfc = (Qik)i for k = 0, • • • , M — 1 and i = 1, • • • , N } and 
the discrete Lagrangian 



M 



£ = EHl V fc- X fcH 2 +[X fc -X fc _ 1 -rF(A fc _ 1 ;a)]-P fc _ 1 

k=i 

+ [Y k - Y fc _! - rG(A*_ i; a)] • Q fc _!} . (17) 

where we have taken H = F + iG. For this Lagrangian, relation (14) is equivalent 
to the direct equation (10), while relation (15) yields the equations 

SC ^ = 0, (18) 



dX ik dY lk 

where these equations hold for all amplitudes which are not fixed by boundary 
or initial conditions, that is for k = 1,---,M and i = l,---,iV. After some 
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algebraic manipulation, conditions (18) may be put in the form 

Pm-i = 2t(Vm — Xm), 

Qm-i = 0, (19) 

(notice cm = 1) and 

P fc _! = P k + r[H 11 P k + H 12 Q k -2e k (\ k -X k )], 

Q fc _! = Q fc + r[# 21 P fc + # 22 Q fc ], (20) 

for < k < M, where we introduce the N X N matrix differential operators 
H n = (dFi/dXj)^, H 12 = (dFi/dYj)^, H 21 = (dGi/dXj)^, H 22 = (dGi/dY^ 
evaluated at the discrete time k. The action of these operators may be written 
explicitly as 

Hn,ij Pj,k = \r — l r (3Xf k + Y? k ) + 2liXi k Yi k ^ Pi k 

+ 9r(Pi-l,k + Pi+l,k — 2Pifc), (21) 

Hu,ij Qj,k = \u — l r 2XikYik — li(3Xf k + Y^)] Qi k 

+ 9i(Qi-l,k + — 2Qik), (22) 

#21, ij = [— ^ — ^r2Xifclijt + /;(3Y^ + ^ 2 fc)] -Pifc 

— 9i(Pi-i,k + -Pi+i,fc ~~ 2Pifc), (23) 
pL22,ijQj,k = \r — l r (3Xf k + Y^) — 2liXi k Yi k ^ Qi k 

+ 9r(Qi-l,k + Qi+l,k — 2Qik)- (24) 

Equations (20) are called the backpropagated equations, since they may be solved 
iteratively by decrementing the discrete time k, using equations (19) as initial 
conditions. The gradient of the Lagrangian with respect to the parameters, which 
by (16) is also the gradient of our cost function, reads 



da * ~ h 



M-l 

/ / r I AA 7. ' >t \ / / t -M-l AA 7. ' >t \ 

(25) 



' aiXA^a) <9G(A fc ;a) 
"k - 7. r 14k ■ 



da; da; 



This gradient may be computed in the following way: (i) For a given value of the 
parameters and initial conditions, compute the A k from the direct equations (10). 
(ii) Compute the and from the backpropagated equations, (iii) Compute 
the gradient using (25). 

It is possible to generalize this procedure to the yield the gradient of the dis- 
tance J with respect to other quantities. For instance, the imaginary parts of 
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the initial amplitudes may be considered as parameters of the problem. Differ- 
entiating the Lagrangian, one obtains 



dY tl 



dC 



T 



1 



Q 



iO ~ 



dGj(A ;BL) 

dY t0 



Q. 



jo - 



dFj(A ;BL) 

dY t0 




(26) 



Likewise, it would be possible to consider parameterized boundary conditions 
to replace (6) although we have not performed this kind of optimization in the 
work reported here. Similarly, we always use X = V for the real part of the 
initial conditions, although the large error in the data acquisition may justify an 
optimization step for X as well. When the optimization algorithm is not used 
to pick Y , we use a Hilbert transform to recover the imaginary part, as in [5]. 
Specifically, for a signal xit) which is assumed to be the real part of an analytical 
signal z x (t) one has [16] 



where V.P. is the Cauchy principal value. This transformation is strictly valid 
when the time series has infinite length. For a finite-length time series obtained 
from simulations, our numerical experiments have shown that the Hilbert trans- 
form poorly converges to the true imaginary part of the amplitudes. However, if 
the Hilbert transform is performed on the whole data series of length 256, but the 
first and last 30 points are excluded, a reasonable approximation to the imaginary 
parts is obtained. It is worthwhile to stress this last point, since we have found 
that it is not possible to obtain satisfactory results without this elimination of a 
string of data at the beginning and the end of the time series. 



We present here some numerical results pertaining both to experimental data and 
to artificial or "synthetic" data. The optimization procedures for model (4) were 
performed in various ways. In a first procedure the initial imaginary parts of the 
amplitudes are obtained by the Hilbert transform, described above, of the real 
amplitudes. We have varied the length of the training window (t 0} t + jm&xT) 
and found that the best reconstruction was attained near jmax = 30. This 
corresponds approximately to 2,5 periods of the basic solution and is of the order 
of the correlation time of the system. On Figure 3 we plot the value of the 
optimal parameters as the time step r of the integration is varied. We observe 
that the parameters have a well-defined limit as r/t — > 0. These parameters for 



z x (t) 



x(t) -lTH[x(t)] 



where 




5 Results 
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Table 1: Estimated values of the parameters 





30 — 60 training window 


average for all windows 


8 parametersaverage 


rT 


-0.0084 


-0.0060 ± 0.003 


-0.019 ± 0.016 


ujT 


0.528 


0.535 ± 0.006 


0.50 ± 0.007 


g r T 


-0.00523 


-0.0045 ± 0.001 


-0.0041 ± 0.0006 




-0.00259 


-0.0038 ± 0.0015 


-0.0043 ± 0.001 


l r T 


0.01707 


0.029 ± 0.01 


0.15 ± 0.16 


UT 


-0.0986 


-0.081 ± 0.02 


-0.37 ± 0.05 


d r T 






0.41 ± 0.4 


d t T 






0.46 ± 0.07 



the first computed window are shown in Table 1. This first window is the 30 — 60 
window because of limitations arising from our use of the Hilbert transform. The 
parameter values are made non-dimensional using the sampling time T = 40ms. 
Notice that the parameter r is negative, but that the amplitude B of plane 
waves in phase opposition (q = ir) in equation (8) is still real. The parameter 
values obtained in the first (30 — 60) training window are used in a numerical 
simulation of the model (4), which yields the reconstructed data. Figure 5a shows 
the reconstructed data (dotted line) and the experimental data (solid line). The 
series is plotted over the length of the 30 — 60 window. To quantify the error we 
define the normalized prediction error e k : 

e k = \\X k -\ k \\ 2 /\\(\ k -(\ k ))\\ 2 

The normalized error e k may be interpreted as follows: when e k = the predic- 
tion equals the data while, when e k = 1 the prediction does not approach the 
data better than the average of the signal. The normalized errors for the first 
training window are small, indicating that our fit is very satisfactory. Figure 4a 
present the entire simulation over the full duration for which we have experi- 
mental data. Figure 4b is the corresponding power spectrum. We also checked 
that the characteristic appearance of oscillatory periods with interspersed low 
amplitude or "laminar" phases is recovered even at times much longer than the 
correlation time (see Fig. 12b). On the other hand, after the initial correlation 
time, the error now grows to (9(1) levels. This growth of the error is unavoidable, 
as it is related to the sensitive dependence to initial conditions. 

Another qualitative aspect is that we recover the same phase relationship as in 
experimental data: one observes that wakes oscillate in phase opposition (q = ir 
in equation (7)). This is in stark contrast with the results of the POD analysis 
of [2] where neither short time predictions nor phase relationships are recovered. 
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We also attempted a second, rather different procedure for parameter esti- 
mation. Instead of estimating the initial imaginary parts using the Hilbert 
transform we included, using expressions (26), these imaginary parts in the set 
of parameters for which we try to optimize. This yields a new estimate of both 
a and Yi . Figures 6 shows the long-time predictions and the short-time error 
after optimization of these 22 parameters. The prediction error over the training 
window (t , t + jm&xT) is somewhat smaller. On the opposite the long time 
simulation (Fig. (6)a) differs drastically from the data: several cylinders have 
much smaller amplitudes than in the experiment. This indicates a likely over-fit 
situation. For the results reported in what follows, we reverted to the estimation 
of Yio using the Hilbert transform. 

It is also interesting to ask whether we can distinguish between several func- 
tional forms for the equation (4). In a preliminary attempt, we extended our 
model by adding a higher-order nonlinear term. We changed H to 

Hi(A; a) = (r + iw)A t - + g(A i+1 + A_ x - 2A 8 ) - Z|A,fA t - - rf|A,-| 4 A t - (27) 

with a new complex parameter d = d r + ic/;. We found that the cost function 
has many relative minima. This makes the optimisation procedure much more 
difficult, As a preliminary result, we selected the data in the window (30,90) in 
which a global minimum is easy to find. This yields the parameters in column 
4 of Table 1 and the predictions shown on Figures 7 and 8 . The short-time 
error is approximately the same. The long-time results show a qualitatively 
different behavior, with much longer intermittent episodes. Thus it seems that 
the simplest, 6 parameter model is better. 

A test of the consistency of our reconstruction is obtained by shifting the 
starting point t of the training window. Figure 9 shows the variation of the 
optimized parameters as a function of t . There is a great deal of scatter in this 
result. Since the reconstructed signals are always satisfactory, we interpreted this 
scatter as a degeneracy in the optimized system. To give a simple example of such 
a degeneracy, if our data were a space-and-time-periodic signal of the form (7), it 
could be reconstructed with parameters in a codimension-2 set, constrained only 
by the 2 relations (8) and (9). In other words, all parameters in a 4-dimensional 
manifold would fit periodic data. This degeneracy can be lifted only by non- 
periodic episodes in the data. Such episodes however are relatively sparse in our 
data set. Moreover there is still the possibility of a lower order degeneracy for 
more complex solutions than the plane waves (7). 

To test our optimization method independently of the experimental data set 
we produced a set of synthetic data. We attempt to mimic as closely as possible 
the experimental process, (i) First we simulate the model equation (4) with the 
parameters of Table 1, column 2. We already know that these parameters yield 
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a behavior similar to that of the data. This yields a series of amplitudes A^, 
from which we extract the real parts Xik- (ii) We then average the signal over 
the sampling period T to approximate the time-averaging produced by the video 
recording, (iii) We then "round-off" or "discretize" the output signal, by letting it 
take only a small number of integer values. This round-off mimics the processing 
of the digitized images, in which the wake positions are located at an integer 
number of pixels away from the centerline. The first two steps result in the signal 
Vik defined for every t°^ s 

V}k ] = jJ2 X ^L+n + W %k (28) 

n=l 

where L = T/t and Wik is a decorrelated noise. The third step is 



V, = E 



S 



v + 0.5 

L nnax 



(29) 



where E(x) is the largest integer smaller than x and Vmax is an upper bound for 
the simulated data Vik- The number S controls how discretized the final signal 
is: the resulting data Vik may have integer values in the interval ( — S, S). This 
results in a "deterministic" loss of information, which is quite different from the 
kind of loss of information that occurs with an added noise source. This should 
be compared with the discreteness of the original data, which were equivalent to 
integer values in the range (—5,5). An example of the resulting data and the 
original experimental data at the same scale are shown on figure 10. 

Since the true values a 4 - of the parameters are known for synthetic data, we 
represent the estimated parameters a 4 - divided by a 4 -. Figure 11 shows the varia- 
tion of parameters with the starting point t of the training window. The noise 
level was adjusted to mimic the experimental data. We let the ratio of the noise 
standard deviation to the data standard deviation be ^[Wy /crfX^] = 0.085 . 
Even with such relatively high noise levels, only a fraction of the variability on 
figure 9 is recovered. The parameters are barely affected by the averaging pro- 
cedure (28), and we have found that most of the error comes from the Hilbert 
transform and the threshold process and not from the noise. Clearly we have 
explained only a part of the variability. It is of course possible that the measure- 
ment process introduces a larger noise amplitude crfWy, but other explanations 
are also possible. For instance, one may imagine a slow drift of the parameter 
values during the experiment, 

Finally, we have attempted to optimize parameters for the "large-time-step 
model" obtained by letting r = T in equations (10). Figures 12a and 13 show the 
entire original data set and its simulation where experimental initial conditions 
are used and imaginary parts are calculated by Hilbert transform. Figure 12b 
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shows a long run of the moduli \A\. This plot demonstrates that intermittency 
persists over long times. Similar long time results (not shown here) were obtained 
in the continuous case. It is quite interesting to note that the accuracy of the 
reconstruction in the large-time-step model is of the same order as that of the 
"converged model" in which r <C T . 

6 Conclusion 

We have introduced a method for the optimization of parameters appearing in 
a system of coupled Stuart-Landau equations. The optimized model appears to 
reproduce the data in a very satisfactory way, making possible both short term 
predictions and a longer term reproduction of the qualitative features of the data, 
such as intermittency. The method thus appears promising for the reconstruction 
of dynamical systems from noisy data. 

On the physics side, it appears that the experimental data may be reproduced 
quite accurately by a continuous-time discrete-space Ginzburg-Landau model. 
This is perhaps not surprising given the previous success of qualitative compar- 
isons between such systems and experimental data. However, it is to be noted 
that the experiments are conducted far from the Hopf bifurcation at Re c , where 
it may be surprising that a theory designed for the vicinity of the bifurcation 
may still hold. Moreover, it is interesting to compare our estimated values of 
the parameters with those found in the literature. We notice that the coupling 
renormalizes the real part r of the growth rate: while r is negative in our recon- 
structions, the plane wave solutions with q = ir still grow. It is also interesting 
to compare the ratio / 8 // r to previous numerical and experimental estimates. For 
arrays of cylinders the POD methods of ref. [2, 3] give / 8 // r = — 1.4 albeit with 
undetermined accuracy. For a single, infinite aspect ratio cylinder near threshold 
(see [17], [18] and references therein), / 8 // r = —2.7 ± 0.1 was found but this ratio 
seems to decrease towards —1 as the aspect ratio decreases. We find / 8 // r = —2.8, 
from column 3 of Table 1 but as for the result of [2, 3] we may have a very large 
error. It is interesting to go byond the fit of parameters for the coupled Lan- 
dau equation to investigate other models. Our results in this direction are only 
preliminary. Higher order nonlinearities as in (27) make optimisation much more 
difficult. They also produce qualitatively different solutions, with longer intermit- 
tent episodes. Recently, P. Legal and his coworkers have observed steady states 
in which the wakes of some cylinders stop oscillating indefinitely. We obtained 
this behavior for some parameter sets of the extended model (27). 

Looking further to the general applicability of our method, It is interesting 
to discuss the dual nature, spatial and temporal, of our data and models. In 
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the method we have used, nothing requires in principle our data to be extended 
spatially. Thus the method could be used just as well for dynamical systems with 
a small number of degrees of freedom. On the other hand, the spatial extension of 
our system results in the production of much more data than in a small system. In 
fact, it is possible in theory to increase the number N of interacting oscillators in 
such systems while keeping the number of parameters constant. In fact, systems 
in higher dimensionality should be more interesting from that point of view. One 
would expect the amount of data to grow exponentially, much faster than the 
number of parameters as the dimensionality is increased. 

One difficulty with the present approach is that the parameters cannot be 
unambiguously determined as they vary with the training window used in the 
simulations. The origin of this variation is unclear, although part of it may be 
recovered in trials performed with synthetic data. Determining the cause of this 
phenomenom is certainly a worthy topic for future research. 

The method we describe in this paper may also be a basis for the more am- 
bitious goal to optimize to find the equations of motion themselves rather than 
simply parameters of a given equation. We have tried to introduce a modest 
variation of the model by adding a higher order nonlinear term, but a much more 
systematic search for better models should be possible. We believe that a more 
thorough experience with such optimization techniques would open the door to 
model optimization for a large number of spatially extended systems. 
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Appendix 



In this appendix, we shortly describe the PLMA algorithm. Let a be a given 
initial guess for the parameters, where the sought parameters minimize the cost 
function J {a). Let k denote the current iteration, starting with k = 0. Each 
iteration requires the gradient vector Gk evaluated at the &th estimate of the 
minimum a^. A vector dk (know as the direction of search) is computed and the 
new estimate a^+i is given by + a^dk where the step length is chosen to 
minimize the function J{a -\-a,k dk). A quasi-Newton method is used to compute 
the search direction dk by updating the inverse of the approximate Hessian (Hk) 
and computing 

dk+i = —Hk+iGk+i (30) 
The updating formula for the approximate inverse is given by 

Hk+i = Hk J— (H k y k sl + SkylHk) + -f—(l + Vk T kVk ) s k s l (31) 

y k s k yi sk y{ sk 

where y k = G k +\ - G k and s k = a k+ i - a k = a k d k 

The algorithm uses a two-step method (with H equal to the identity matrix) 
described in detail in [15] in which restarts and pre-conditioning are incorporated. 
The termination criterion involves a pre-defined tolerance e T . The algorithm is 
halted if the following three conditions are satisfied 

i Jk-l-Jk < Cr(l + \Jk\) 

ii ||a fc _i - a k \\ < e l J 2 (l + ||a fc ||) 

iii HGfcll < ey 3 (l + \ Jk\) or 1 1 1 1 < e^, where ca is the absolute error associated 

with computing the cost function J . 
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Figure 1: a) Snapshot of the experiment, b) Original data acquisition for fixed 
time. Peaks locate tracer lines. The deviation of tracer lines from their baseline 
position defines the experimental wake amplitudes Vik- 
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Figure 2: (a) Experimental data composed by 256 original acquisitions. (Showed 
in the range 30-256) (b) Power spectrum of the original signal computed over the 
identification window. The inset shows the spectrum computed over the entire 
time series. 
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Figure 3: The parameters r and l r defined in the text as identified with various 
values of the time step r. The smallest integration time is 2 10~ 4 s which gives 
2400 integration points by oscillation period. 
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Figure 4: (a) Simulation of the continuous model. The parameters used cor- 
respond to the smallest r in Figure 3. Initial conditions are provided by the 
experimental data, (b) Power spectrum of the simulated signal computed over 
the identification window. The inset shows the spectrum computed over the 
entire time series. 




Figure 5: (a)Over plot of the coupled Landau equation simulation (dotted line) 
and experimental data (solid line) (b) the normalized error as defined in 

the text. 
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Figure 6: (a) Simulation over the full range of the data set using the parameters of 
the 22 parameters optimization, (b) Short time normalized error in the training 
window. 




Figure 7: (a) Over plot of the simulation with the extra term |A| 4 A added to 
the coupled Landau equation in the hrst training window. The dotted line are 
simulations and solid lines are experimental data. The parameters averaged from 
training windows from 30 to 90 are used, (b) Short time normalized error in the 
training window. 
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Figure 8: Simulation over the entire duration of the data set with the extra term 
\A\ 4 A. 
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Figure 9: Parameters r and l r as function of the starting point of the optimization 
window. The identification is performed in a window of 30 time steps taken from 
experimental data set. 



vAyvwwwwwvA 




vA\y\yv\A7V\AAAAAAAAAAAAAy 
vwwvWvA/VWVVWvW 



23 



ORIGINAL DATA 





Figure 10: (a) Experimental data showing the discreteness of the signal due to the 
measurement process, (b) Synthetic data where the discretisation or "round-off" 
parameter S is set to 6. 
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Figure 11: Normalized parameters r ■normalized = r/r and ^normalized = 
l r /l r as functions of the starting point of the optimization window when the 
round-off parameter S = 6 
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Figure 12: (a) Simulation of the large-time-step model where the integration 
timestep coincides with the measurements (r = T) and (b) The moduli \A{\ in a 
simulation pursued over times much longer than the duration of the experimental 
series. Note that the intermittency, characterized by holes in the moduli persists 
over the entire simulation. 
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Figure 13: (a) Over plot simulation (dotted line) for the large-time-step model 
and original data (solid line) (b) Normalized error for the optimization window. 
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